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Abstract 

A Navier-Stokes procedure to calculate the flow about 
an airfoil at incidence has been developed: The parabolized 
equations are solved in the streamline coordinates generated 
for an arbitrary airfoil shape using conformal mapping; A 

modified turbulence model is applied in the entire 
domain, but the eddy viscosity in the laminar region is 
suppressed artificially to simulate the region correctly. The 
procedure has been applied to airfoils at various angles of 
attack and the results are quite satisfactory for both laminar 
and turbulent flows. It is shown that the present choice of 
the coordinate system reduces the error due to numerical 
diffusion and that the lift is accurately predicted for a wide 
range of incidence. 


Introduction 

Aerodynamic characteristics of an airfoil at incidence, 
especially near and beyond the stall angle, is of paramount 
practical interest as these are closely related to the 
performance of engineering devices such as aircraft and 
turbomachinery. Because of the importance associated with 
the flow, much efforts have been devoted to develop 
prediction techniques for these flows. 

One may approach the problem by using the interactive 
methods that explicitly couple the viscous and inviscid 
effects in an iterative manner. The methods by Maskew & 
Dvorak, 1 Gilmer & Bristow, 2 and Cebeci et al. may belong 
to this category. By and large, the methods have been 
successful in predicting C JmMX and the subsequent stall. 
However, since these all adopt the boundary-layer 
procedure, special treatments are necessary to handle the 
reversed flow region; the details of the flow or the accuracy 
in this region may suffer. 

On the other hand, the method based on the Navier- 
Stokes equations, which is gaining popularity with the 
advent of modem computer technology, is more rigorous 
and appropriate, in principle, than the former for the highly 
interacting flows as the equations are valid both in potendal- 
and viscous-flow regions. Handling of the separated region 
is more straightforward, too. Among many earlier attempts, 
Shamroth & Gibeiing 4 made compressible-flow calculations 
using a transitional k< turbulence model and Rhie & Chow 5 


used a SIMPLE-type method with a k-e model to predict 
the pressure distribution and the near wake. Both adopted 
nonorthogona! computational grids due to their simplicity 
and generality. Although these calculations exhibit certain 
degree of success, the results are not entirely satisfactory: 
the flow at near- or post-stall angle has not been 
successfully predicted. Chang et al. 6 observed the similar 
shortfalls in an existing Navier-Stokes procedure in their 
comparative study of interactive boundary-layer and thin- 
layer Navier-Stokes procedures. 

The intention of this paper is to present a new Navier- 
Stokes procedure, in which the various aspects of the 
calculation have been improved, and to show that the flow 
over a wide range of incidence has been predicted with 
reasonable accuracy and robustness. 


Grid Generation 

Among various grid generation techniques, a method 
based on conformal mapping has been adopted as it has 
distinct advantages in treating the flow of present interest. 
Specifically, the grid lines so generated are orthogonal to 
each other and, moreover, the coordinates can readily be 
made to be intrinsic. These two points are not imperative. 
However, the equations do become simpler when the 
coordinates are orthogonal and the false diffusion in the 
numerical scheme is greatly reduced if the coordinate line is 
aligned with the local streamline. The results are, therefore, 
expected to be more accurate. 

The conformal mapping used here transforms an 
arbitrary airfoil shape onto a unit circle by two successive 
transformations. The profile in the physical plane z is first 
transformed into a smooth near circular section in the plane 

£ by the Karman-Trefftz transformation, which removes the 
sharp comer at the trailing edge, and, subsequently, into a 
unit circle by solving the uershgorin integral equation. The 
latter part is done numerically after the integrand is suitably 
modified to make the procedure more tractable and accurate. 
The details are referred to Choi & Landweber and will not 
be repeated here. 

The resulting mapping relations may be written as 

z=flO 0) 

and the Laurent series. 
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where/ denotes the Karman-Trefftz transformation and A & 
a's are coefficients that are determined from the second 
transformation. It is important to point out that, since the 
profile in the intermediate plane £ is nearly circular, the 
number of terms required in the series, Eq. (2), to accurately 
compute C is not large: 10 terms have been found sufficient 
and used in the present work. 


because of nonzero circulation, T, the potential at the trailing 
edge is double valued. 

<fru = < t > rL + r ^ 

and a jump in <p is present across the trailing streamline. 


From these relations, various types of grid , i.e., C-, 
H- and O-grid, can now be constructed. The radial lines 

and the concentric circles in the t plane give an O-type grid 
while the horizontal and vertical lines in the plane of 
complex potential W) and those in the VV plane give, 
respectively, H- and C-type grids. The grid of H-type is 
used in the present calculation and the details of how it is 
obtained is described below. 

The complex potential W for a stream velocity U at an 
angle of attack a about a unit circle at the origin is 

tv = ” J +i (3) 



where T is the circulation about the circle and is equal to 
4nUsin(a-d°) so that Eq. (3) satisfies the Kutta condition 
that the velocity be zero at the trailing edge, B—0 o . The 
velocity U in the T plane is related to the undisturbed 
velocity U„ in the z plane by 
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Fig. 1 Physical and computational domains for the flow 
about an airfoil at incidence. 


The coordinate lines in the W plane are lines of constant 
potential (<p) and stream function (iff)', corresponding lines in 
the physical plane are also equipotential lines and 
streamlines of the flow under consideration, and constitute 
L orthogonal grid of H-type. One point to observe is that. 


(b) 

Fig. 2 Sample grid and the close-up view of the nose region 
about an airfoil at a=5 deg. 


The remaining task to complete the grid construction is 
:o distribute the grids efficiently. This is accomplished by 
using tank as a distribution function to place more grids 
where needed, e.g., near the surface, around the leading and 
trailing edges. For the proper clustering in the streamwise 
direction, the grids are First distributed along the stagnation 
streamline ABCD and AB'C'D shown in Fig. 1 using the 
arc length as parameter. The number of grids for the 
segment BC may be different from that for B 'C . The 
transformed grids in the W plane can then be obtained by 
using the relations (1), (2), and (3). However, a direct 
attempt to do so involves rather time-consuming algebra; 
the following spline interpolation is used instead. For a 
given set of points along the <p axis, the corresponding 


points, which lie on the stagnation streamline, in the z plane 
are readily determined: r from Eq. (3) by Newton’s 

rootfinding algorithm, f by Eq. (2) and z by Eq. (1). The 
arc length, $, for each of these points is then calculated and 
the relation between s and <p is established. A cubic spline 
function is used to relate the two and, for a point in the z 
plane, this interpolation function gives the matching point on 

the <t> axis in the W plane. The grid clustering in the vertical 
direction, on the other hand, is done in the W plane using Y 
as parameter. A typical grid in the physical plane for ct “5° 
is shown in Fig. 2. 

Governing Equations 

Following Nash & Patel , 8 the continuity and Reynolds- 
averaged Navier-Stokes equations in general orthogonal 

curvilinear coordinates (£,tj ) are written as: 

Continuity: 

T^{-k ih 2 u ) + 'k {h ‘V)} =0 (6) 
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where (U f V) and (u,v) are the mean and fluctuating velocity 
components, respectively, in the (£t) direction, p the 

Uj: 

pressure. Re ( =— ) the Reynolds number, v the 

kinematic viscosity, and h and K the metric coefficients and 
curvature parameters. The equations have been made 
dimensionless by using the freestream velocity and the 
airfoil chord c . These equations are of conservative form 
and are exact except for the neglected streamwise diffusion 
terms. The conservative form appears to give more stable 
behavior of the numerical method in the neighborhood of 
the stagnation point where the H-grid becomes singular. 

The Reynolds stresses in Eqs. (7) and ( 8 ) are related to 
the mean rates of strain through the eddy-viscosity 
hypothesis and are given in the next section. 
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Turbulence Model 

A modified k-e model is adopted as a closure 
relationship in the present study. The transport equations 
for the turbulent kinetic energy and the rate of dissipation 
(7) compatible with Eqs. ( 6 )-( 8 ) are 
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Solution Procedure 


and the model constants C u . <7 k , a c , C £l , C e2 andC eJ are 
given the values of 0.09, 1.0, 1.3, 1.44, 1.92 and 4.44, 
respectively. 

The dissipation equation was first proposed by Hanjalic 
& Launder, 9 where they reasoned that the energy transfer 
rates across the spectrum are preferentially promoted by 
irrotational deformations and showed the improved 
prediction over the standard k-£ model, especially in the 
adverse pressure-gradient region. It should be noted that 
the dissipation equation assumes the | direction to be the 
predominant flow direction and the present intrinsic 
coordinate system is consistent with this assumption. 

Two key modifications to this model have been made 
for the present study. Rather than using the wall function in 
the near-wall region, the two-layer approach of Chen & 
Patel 10 is adopted to make the model applicable in the 
separated-flow region and to provide a finer resolution in the 
near- wake region. The other change made is the use of the 

anisotropic k-£ model of Nisizima & Yoshizawa 11 to 
represent the Reynolds normal stresses. The Reynolds 
stresses are then expressed as 

-™= v '{t,% + t 2 %- k ‘ 2U - K2iV ) (13) 

= (14) 

-v 2 = ~jk + 2 v, + K 2I U ) + S v i (15) 

where 

S.,,L(-2C„+C«)£{-L %-w) 2 (16) 

v4<c„-Jc a )L(_L %-w) 1 (17) 

and C T/ = 0.07, C x2 - -0.015 from Ref. 11. The 
nonlinear terms in Eqs. (14) and (15) lead to the anisotropy 
of the turbulence intensities. 

When using the k-e model, it is customary to assume, 
for computational convenience, that the flpw is turbulent 
everywhere as was done in Rhie & Chow. Although this 
may be justifiable as the laminar portion of the flow is 
limited to the small region near the nose, the artificially 
produced turbulent flow, which may be healthier than the 
real flow, could greatly affect the leading-edge-separation 
pattern. 

In order to get around this difficulty inherent to the k-e 
model, the concept of intermittency is employed: the 
transport equations for k and e are solved in the entire 
domain, but the eddy viscosity is set to be zero in the 
laminar region. The procedure has been found more 
successful than solving the equations only in the turbulent 
region. The latter performed relatively poorly as the initial 

profiles for k and £ at the transition location could not be 
provided accurately. 


The governing equations, Eqs. (6)-(10), are solved in 
the calculation domain bounded by constant £ and Tj lines. 
Using the staggered grid, shown in Fig. 3, the diffusive 
derivatives of the equations are discretized by central 
differencing while the convective derivatives in the 
streamwise and cross-streamwise directions by upwind and 
hybrid differencings, respectively. The numerical scheme 
adopted in the study is the modified version of the CELS 
(Coupled Equation Line Solver) algorithm used in Ref. 12: 



Fig. 3 Grid layout and storage location for each variable. 


The calculation proceeds in the streamwise direction 
and the solution at a given streamwise station is obtained 
simultaneously. The penta-diagonal system of equations for 
V, which is derived from the continuity and_ momentum 
equations by eliminating the pressure and the streamwise 
velocity component, is solved first. The pressure and U then 
follow successively by the backward substitution. Using 
these values, the turbulence transport equations are solved 

for k and e by the Thomas tridiagonal matrix algorithm. To 
enhance the convergence, a backward pressure correction is 
applied at the end of each complete sweep. This is 

accomplished by forcing the | -momentum equation be 

satisfied, on the average, along ea ch co nstant £ line. The 
process is repeated until the specified convergence criterion 
is met. Maximum pressure variation of 10 is used for the 
present calculation. 


The calculation is performed for a sufficiently large 
domain that encompasses the entire profile and the following 
conditions are specified at the boundary: 

dV _ 

upstream: U = Up# * ~ 0 


downstream: 

outer 



U-U^. 


dV 


= 0 


wall: 


no-slip condition 


where the subscript^/ indicates the potential-flow value. 
The turbulence quantities at far boundaries except along the 
downstream end, where the conditions on turbulence are not 
needed, are assinged a very small value to simulate the non- 
turbulent flow. 


Results and Disscussion 

Laminar flow 

The calculation is first performed for the laminar flow 
about a 12%-thick symmetric Joukouski airfoil section at 
/?e=1000. For the incidence angle of 5°, the grid of 
(140x40) and the calculation domain which covers the 
region -1 < x/c < 5, -3 < y/c < 3 were found adequate. A 
coarser grid (70x40) appears to give comparable results and 
an optimum grid may lie somewhere in between. However, 
no further attempt has been made to find this grid 
distribution. 

The velocity vectors and the surface pressure 
distribution are presented in Figs. 4 and 5. The flow 
separates at about midsection and, consequently, the 
pressure distribution is altered substantially from that of the 
inviscid flow. The results are seen to be in excellent 
agreement with those by Ghia et al. 13 , who solved the 
streamfunction-vorticity equations on a (229x45) C-type 
grid. 

Figure 6 illustrates the importance of the grid alignment 
with the flow. Here, the calculations have been made for a 
= 8° with two different grids: one grid is generated for a - 

0° and the other for a = 8° and, as a result, the former is 
skewed by 8° in relation to the flow direction. It is observed 
from the figure that a finer grid is required when the grid is 
skewed to obtain the results of comparable accuracy. It is 
primarily due to the numerical diffusion caused by the first 
order upwind differencing and the discrepancy could be 
reduced by incorporating a higher order upwind scheme. 
This will, however, introduce additional complexities into 
the coding and it is desirable to construct a grid which 
follows the general flow direction whenever possible. 



Fig. 4 Velocity vectors for the 12%-thick Joukowski airfoil 
section at a=5 deg and Re =1000. 



Fig. 5 Surface pressure distribution on the 12%-thick 
Joukowski airfoil section for a =5 deg and Re- 1000. 



Fig. 6 Comparison of two different grids for a=8 deg. 




Turbulent flow 

For turbulent flows, the calculations have been 
performed for NACA airfoil sections, namely 4412 and 
0012, at various angles of attack. A 140x40 grid is fitted 
over -1.5 < x/c < 10, with the first point normal to the 
surface being placed approximately at y + *5 . It is 
reminded that the grid needs to be reconstructed when the 
angle of attack or the Reynolds number varies. The vertical 
boundary is located at about where the tunnel wall is to 
closely mimic the experimental condition and the slip 
condition is imposed there. Since the wall and the constant 

7} line do not coincide, we introduced a vertical velocity 
component of right amount during the computation to make 
the velocity vector parallel to the tunnel wall. The wall 
location is indicated by the dotted line on the present grid for 

the NACA 4412 airfoil at a = 13.9° in Fig. 7. 



Fig. 7 Computational grid with tunnel wall location for the 
NACA 4412 airfoil at a=13.9 deg. 


The pressure distribution for the NACA4412 airfoil at 
a - 13.9° and Re = 1.5xl0 6 is presented in Fig. 8: In the 
calculation, the transition points for upper and lower 
surfaces are prescribed to be at 0.025c and 0.103c, 
respectively, as done in the experiment. The present result 
(solid line) is in near exact agreement with the experiment of 
Coles & Wadcock. As in the laminar case, the pressure 
distribution is altered greatly from the inviscid one by the 
flow separation. Also shown in the figure are the results by 

the standard k- e model with and without the tunnel wall 
effects. Here, the result without the tunnel wall means that 
the calculation is performed in a larger domain (-5 < ylc < 5) 
with the freestream boundary condition. Although these all 
are in relatively good agreement, it is evident that each of the 
changes results in nodcible discrepancies. 

Figures 9 and 10 show the velocity vectors and the wall 
shear- stress distribution. The velocity vectors and the grid 
lines, which are the streamlines of the inviscid flow, 
coincide in most of the region. This is expected and 
validates the present choice of turbulence model and the 
approach of parabolization. The wall-shear stress shows 
that the boundary layer separates at xlc *■ 0.8 and the 
laminar boundary layer is very close to separation before it 
becomes turbulent at xlc = 0.025 . It is cautioned here, 
however, that the absolute values of wall-shear stress in the 
upstream section of the airfoil may not be accurate as the 
boundary layer is too thin to be adequately resolved by the 
present grid distribution. To check how well the turbulence 
model mimics the transition process, the turbulence 
quantities, k and v r , and the wall- shear stress in the 
neighborhood of the transition point are examined in Fig. 
11. The turbulent kinetic energy and the eddy viscosity 
plotted are the maximum values at the given station. The 
smooth but rather sharp increase in these quantities indicates 
that the present treatment for transition is qualitatively 
correct The turbulent kinetic energy does not grow in the 
laminar region because the prodution terms in the transport 
equation are turned off by suppressing the eddy viscosity. 




Fig. 9 Velocity vector for the NACA 4412 airfoil section 
fora=13.9 deg and/?e=l. 5x10 . 


Fig. 8 Pressure distribution on the NACA 4412 airfoil 
section for ct=13.9 deg and Re =1.5x10 . o, 

experiment 1 4]; , present; , standard 

jfc-£ model ; — — , standard k ~€ model with 
freestream condition; — , inviscid flow. 






Fig. 13 Pressure distribution on the NACA 0012 airfoil 
section for oc=6 deg and Re = 1.5x10 : O , 

experiment[15]; , present. 



Fig. 14 Pressure distribution on the NACA 001 £ airfoil 
section for a =6 deg and Re = 2.8x10 : O , 

experiment[15]; , present ; * Shamroth 

& Gibeling[4]; , Rhie & Chow[5]. 


TheC, - a curve for/?e = 1.5xl0 6 is depicted in Fig. 

15, along with the curves obtained by the standard k-e 
model with and without the wall effect. The results are in 
good agreement with the data and show the similar 
characteristics as in the case for the NACA 4412 airfoil. 
For the incidence angle greater than that shown in the figure, 
the anisotropic turbulence model becomes less stable due to 
its nonlinear terms; the convergence is slowed as more 
under-relaxation is required. The calculation was thus not 
carried out for the case much beyond the stall angle. 



Fig. 15 C r a curve for the NACA 0012 airfoil section at 

Re- 1.5xl0 6 : A, experiment 15] ; , present ; 

, standard k-e model ; , standard k -e 

model with freestream condition. 


Concluding Remarks 

A new and improved Navier-Stokes procedure has been 
developed and applied successfully to the flow about the 
airfoil at incidence: the lift of the airfoil is accurately 
predicted for a wide range of angles of attack. It has been 
shown that the present choice of the coordinates, i.e., the 
streamlines and the equi-potential lines of the inviscid flow, 
helps make the method more accurate and efficient. The 

modified k-e turbulence model, which is used in the whole 
domain with zero inteimittency in the laminar region, gives a 
qualitatively correct transition behavior. 
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